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SUMMARY 

A numerical solution to the Navier-Stokes equations is obtained for 
blunt axisymmetric entry bodies of arbitrary shape in supersonic flow. These 
equations are solved on a finite-difference mesh obtained from a simple 
numerical technique which generates orthogonal coordinates between arbitrary 
boundaries. The governing equations are solved in time-dependent form using 
Stetter's improved stability three-step predictor-corrector method. For the 
present application, the metric coefficients were obtained numerically using 

fourth-order-accurate, finite-difference relations and proved to b e to tally 

reliable for the highly stretched mesh used to resolve the thin viscous 
boundary layer. Solutions are obtained for a range of blunt-body nose shapes 
including concavities. Results indicate that the numerically generated 
coordinate system performed exceptionally well and no problems were encountered 
in the coupling of the numerical coordinate generator and the fluid dynamic 
equations . 


INTRODUCTION 

One of the major problems retarding the rapid development of computa- 
tional fluid dynamics for complex geometries has been the difficulty of 
generating the finite-difference mesh. Much effort has been expended to 
develop coordinate transformations and/or mesh generators for varying degrees 
of geometric complexity (see refs. 1-5 for representative examples). For 
viscous flow over blunted bodies, such as planetary entry vehicles, the 
boundary conditions on the body and the large gradients adjacent to the 
surface must be represented accurately by the finite-difference approximations 
to the Navier-Stokes equations. Toward this goal, almost all numerical solu- 
tions to the Navier-Stokes equations to date for blunt-body flows have used 
body geometries conducive to use with natural or nearly orthogonal coordinate 
systems. 

In the natural coordinate system, the body surface itself forms one 
boundary, i.e. the body contour coincides with a constant coordinate line. 
Typical examples of this approach are a cylindrical coordinate system to 
describe flow over a cylinder, a spherical coordinate system to describe flow 
over a sphere, and a parabolic coordinate system to describe the flow over a 
paraboloidal body. Reference 6 gives a representative use of a natural 
coordinate system for the numerical solution of a fluid flow problem. In the 
natural coordinate system, the normal coordinates intersect the body orthog- 
onally, thus simplifying the boundary conditions. There is little difficulty 
in compressing the mesh near the body because the computational mesh system 
is composed of lines parallel to the body which can be concentrated as close 





to the body as desired. There is, however, one rather severe restriction on 
the natural coordinate system; that is, the body must have an analytic shape. 
Unfortunately, most planetary entry vehicles bear little resemblance to the 
limited number of natural coordinate systems available. 

Another option is similar to the natural coordinate system in that the 
body surface becomes one coordinate line in the system. This is called the 
body-oriented coordinate system. In this system, the coordinates of a point 
are determined by the distance along a body surface measured from the axis of 
symmetry and the distance along a normal to the body. This type of system 
has been used to describe the flow over the forebody portions of blunt entry 
bodies; references 7-9 give some representative examples. Although this 
system has received wide use, it does suffer from several major deficiencies, 
one of the most serious of which is that it cannot handle body shapes with 
concavities (ref. 10) . 

Both conformal (ref. 5) and near-conformal mapping (ref. 4) have been 
used to generate coordinate meshes about complicated geometric shapes; however, 
these techniques are mathematically complicated and generally require multiple 
transformation steps leading to a loss of physical reality in the computational 
plane. Such complications make finite-difference mesh setups difficult and 
the computer codes generally are not easily applied/converted to general shapes. 

Recently, reference 11 presented an application of a simple numerical 
coordinate generator (ref. 12) to blunt-body shapes. This technique allows 
for the numerical generation of general orthogonal coordinate systems about a 
wide range of body geometries. In this technique, the body can be represented 
by a series of discrete (but continuous) points rather than by an analytical 
approximation, and in the transformed computational plane, the region of 
interest is rectangular with the body surface being a coordinate line. This 
representation combines the advantages of the natural coordinate system and 
the body-oriented coordinate system. An additional advantage is that the 
coordinates are generated in the physical plane, which simplifies the finite- 
difference mesh setup. This technique was used in a time -dependent solution 
procedure for inviscid flow over blunt bodies (ref. 13) . The use of the time- 
dependent solution procedure, where the moving shock wave was the outer- 
coordinate boundary, required that the coordinate system be regenerated at the 
end of each time step. This coupling of the fluid dynamics and the coordinate 
system worked very well for a large number of test cases. However, this 
application was for inviscid flow only and the finite-difference meshes were 
nearly equally spaced. 

The present analysis uses the coordinate generation technique (ref. 11) 
along with the time- dependent solution concepts (ref. 13) to obtain solutions 
to the full Navier-Stokes equations for viscous compressible flow over blunt 
bodies. The emphasis is on obtaining fluid-flow solutions for blunt bodies 
of varying nose shapes (complexity) where the mesh must be highly compressed 
to resolve the viscous boundary layer. 
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SYMBOLS 
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J 

K c 

M 

N 


N. 


Pr 
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R c 
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forebody shape coefficients 

Sutherland's viscosity coefficients 

specific heat 

average error 

arbitrary function 

arbitrary vector function 

static enthalpy 

metrical coefficients 

transformed coordinate metrics 

£ direction index 

total intervals in £ direction 

D direction index 

total intervals in ri direction 

spacing parameter 

Mach number 

unequal spacing coordinate parameter 
Prandtl number 
Reynolds number 
pressure 

- - -3 

nondimensional heat transfer, q = q/p V 

J 1 1 ~CO 00 

nose radius, m 

local distance between body surface and outer boundary 
radius of body surface 

cotangent of body angle at tangency point 
nondimensional time, t = t V /R., 

’ 03' N 


T 

U 

u 

V 
S V 

V 

X»Y 


nondimensional temperature, T = T C /V 

Poo °° 

arbitrary vector 

nondimensional tangential velocity, u = u/V ( 
free-stream velocity, m/s 
nondimensional normal velocity, S V = S V/V oo 
nondimensional normal velocity, V = V/V^ 
Cartesian coordinates 


S X, S Y 


Z 

a 

3 

Y 

e 

n 

e„ 


u 

c 

p 

p 

a 


shock location in Cartesian coordinates 
axial coordinate 

arbitrary parameter in Stetter’s method 
shock angle 

ratio of specific heats 
convergence criteria 
transformed normal coordinate 
body-surface angle (see fig. 1) 
local angle (see fig. 1) 
cone angle 

nondimensional viscosity, u = u/u^ 
transformed tangential coordinate 
nondimensional density, p = p/p^ 
radial distance (see fig. 5) 
damping coefficient 


Superscripts : 


dimensional values 


Subscripts : 

oo free-stream quantities 
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METHOD OF ANALYSIS 


Coordinate System 

The numerically generated orthogonal coordinates will be determined from 
the original X,Y coordinate system's description of the body surface and 
shock wave. Taking the origin of the X,Y system as lying inside the body 
to be described, the surface distance 5, which forms one of the transformed 
orthogonal coordinates, can be calculated by defining £ as zero on the -X 
axis (fig. 1) and increasing to unity at the end of the forebody surface. 

Thus 5 is given by 



where r = (X 2 + Y 2 ) 1/2 and 0. = cos' 1 (-X/r ) . As in the analysis of 

reference 13, the shock wave is taken as the outer boundary of the transformed 
coordinate system. On the outer boundary, n = 1 while on the body surface, 
n = 0. The level lines between the outer boundary, shock wave, and the body 
surface can be constructed along straight lines connecting corresponding 
points on the body and shock. Note that the mesh points on the outer boundary 
are not the final mesh points, but initial values used only to set up the 
level lines. The actual mesh points will result from the numerical generation 
of the orthogonal normal lines. The spacing of the level lines is arbitrary, 
however, for viscous flows, the boundary layer must be resolved and the unequal 
spacing relation of reference 14 can be easily applied, 

k c V an - 1 

n j = ” l/AN _ j 
c 


where N. = (j-l)AN and Nj = 1., with AN = 1/(J-1) and K c being the spacing 

parameter (generally 1 < K c < 2) . The larger the spacing parameter K c , the 

more unequal the spacing. Using 'the unequal spacing relation, the level lines 
between all corresponding points on the body and shock can be calculated. The 
relationships for the corresponding values of X,Y can be obtained (see fig. 1 
for geometrical schematic) from 
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where 


-1 


e i = si " IP'i.n-i - Y i,n-o’ /r i,iJ- 


Figure 2 shows the level lines constructed in this manner for a spherically 
capped conical body with K c = 1.01, which gives nearly equal spacing. 

Once the level lines have been determined, the normal lines are 
constructed numerically so that an orthogonal system is defined. The approach 
to the construction of the normal lines is the one given in reference 12 
which uses a simple "predictor-corrector" technique analogous to the 
trapezoidal integration method of numerical integration. In this technique, 
the solution is first predicted from the level line at a known point by using 
the Euler method. Once the predicted point on the next level line is 
obtained, the slope at that point is calculated and a new predicted point is 
obtained using this slope. The actual solution is then a combination of these 
two solutions, i.e. the final X,Y values are an average of the predicted 
and corrected ones. This procedure is illustrated in figure 3. Starting on 
the body, the solution proceeds point by point along a level line until all 
normals on that level have been constructed. Then the solution proceeds to 
the next level and the process is continued until the outer boundary shock 
is reached. Figure 4 shows a typical orthogonal coordinate mesh constructed 
about a spherically capped conical body. 


Coordinate Metrics 

Once the coordinate system is constructed, then the X,Y location of all 
mesh points is known and using the coordinate system depicted in figure 5, the 
metric coefficients can be determined using the nomenclature of reference 15. 

1 ~ i 

q X = p cos q) 

2 ~ 

£ x = p sin <J> 

. 3 

(j> x = z 


z = z(£,r|). The metric coefficients are obtained 
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For an orthogonal system, the metric coefficients h^ ^ , h^ h^ h^ ^ , 

h„ _, and h, - all have to be zero, leaving only the three familiar coefficients 
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When the derivatives in the metric coefficient relation are taken, the following 
metric coefficient expressions are obtained: 


h _ f9p,2 . (9z,2 

h i,i - W + w 

v , 3p 9p %z_ 3z_ 

h l,2 ' h 2,l “ 3n 35 3n 95 


h l,3 h 3, 1 0 

u _ (9p-v2 /- 9z- > 2 

h 2,2 " + ( 95° 


h 2,3 h 3,2 ° 


h 3,3 P 


With the choice of the present coordinate system, two of the three necessary 
orthogonality conditions are identically satisfied, leaving only the h^ 2 

coefficient, which was shown in reference 11 to be negligible for the present 
application. Thus, only the three necessary metrics are left: 
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Since the computational plane (5,n) is an equally spaced rectangular region 
(see fig. 6), the derivatives in the above metric coefficient relations can 
be evaluated using equally spaced, central finite differences. For the 
present analysis, fourth-order-accurate relations (ref. 16) are used in place 
of simpler second-order-accurate finite differences in order to produce smoothly 
varying metric coefficients adjacent to the stagnation line. These fourth-order- 
accurate metric coefficients proved to be totally satisfactory for the present 
analysis even in the regions of high-mesh stretching/compression. 


Governing Equations 

The equations used in the present analysis are the full unsteady Navier- 
Stokes equations in general orthogonal coordinates for laminar viscous flew 
as given by reference 17. These equations in nondimensional form are. 
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The axis of symmetry, stagnation line, is a singular line in the computational 
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domain because h 3 = 0 and u = 0 along the axis. The limiting forms are: 
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V • v = =- 
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In addition to the above conservation equations, an equation of state and a 
viscosity law must be specified. For a perfect gas, the equation of state 
can be given by 

Y-l 


P = 


Y 


ph 


and Sutherland's viscosity equation is used in the form 


if 3/2 

u = c i c 2 + T 

to determine the molecular (laminar) viscosity. 


Boundary Conditions 

Along the body surface for the present analysis, the following conditions 
are imposed: 

u = 0 (no slip) 

v = 0 (no blowing) 

h = h w (constant wall temperature) 



Because the flow analysis is in the unsteady form, both the flow properties and 
shock velocity have to be specified at the upper boundary. This is' accomplished 
using a modified form of the Rankine-Hugoniot relations, a complete development 
of which can be found in reference 13. In order to apply these relations, the 
pressure immediately behind the shock has to be calculated using only 
quantities interior to the computational region. This calculation is 
performed using the continuity equation, energy equation, and equation of 
state with the appropriate derivatives in these equations being expressed by 
backward finite differences. Once the shock velocity is obtained at each 
mesh point along the shock, the shock position is determined from the 
solution of the following equations: 


~ = Sv - cos (tt/2 - 3) 
dt i 


d S Y. 

x 

dt 


S V ± sin (tt/2 - 6) 
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The flow conditions along the outflow boundary are obtained using the follow- 
ing second-order-accurate extrapolation relation for equally spaced data 
(ref. 7): 


f. 

1 


( 2 °fi_l - 6f._ 2 - 4f._ 3 + f._ 4 )/ll + OCAC 2 ) 


Numerical Method of Solution 


Because of the necessity of recalculating the orthogonal coordinate 
system at each time step (due to the shock movement) , it is imperative that 
the overall solution converge as rapidly as possible to minimize the 
additional work. For this reason, Stetter's method (ref. 18) is chosen as 
the method of solution, for it has been demonstrated (refs. 19 and 20) to 
have an expanded region of stability which allows for a more rapid marching 
to steady state. In Stetter's three-step predictor-corrector method, the 
derivatives in the governing equations are expressed as second-order finite 
differences so that the overall solution is of second-order accuracy in the 
space dimensions. Since the steady-state solution is the desired end 
product, the dp/dt term in the energy equation is dropped (refs. 7 and 9) 
and thus temporal accuracy is not maintained. 


In order to facilitate the use of Stetter's method, the governing 
equations are recast in the following simplified form: 


U = 



BU 

3t 


F 


and F represents the righthand side of the 
appropriate equation. 


The three-step process becomes: 
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(step 3) 

(parameterization) 
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After the parameterization step, a damping (or smoothing) of. the solution 
is performed to prevent oscillations in the flow from destroying the solution 
process. This damping is the fourth-order procedure given in reference 21 
and used quite successfully in references 9 and 10. 

Because of the large variation in grid spacing used in the present 
analysis, the local maximum time step based on the CFL condition could not 
be used as was done in reference 20. Instead, a local minimum time step was 
found along each E, = constant line and then the time step was defined as 


i=constant 

This time step proved entirely satisfactory for the range of body shapes 
given in the present analysis. 

The orthogonal coordinate system was numerically generated at. the end 
of each full time step while at the end of steps 1 and 2, only the shock 
position and related metrics were calculated. This is the same procedure 
used in reference 13 where a wide variety of inviscid blunt-body flow fields 
were calculated. Note that in the present study, no attempt was made to optimize 
the fluid dynamic/orthogonal coordinate calculation procedure. 

The solution procedure was assumed to be converged (at steady state) 
when the average error was less than the convergence criteria. The average 
error is defined as 



1-1 J-l 


E = 


i=l j=2 


N+l N 

p. . - p. . 

1.3 1.3 


N 

3 . . 

i>3 


(I-DCJ-2) 

Thus, the solution is converged when 


E < e 

where e is a small value on the order of 10 - ^. 


RESULTS AND DISCUSSION 

The present analysis was used to calculate the flow over a 45° hyperboloid 
using the test conditions given in table 1. The computational grid was a 
15 * 31 mesh with a spacing parameter K c = 1.2 in the u direction, which 
gives a highly compressed mesh adjacent to the body surface. The spacing 
parameter for the E, direction was K c = 1.15, which gave a slight 
compression in the stagnation region and a stretching of the mesh toward the 
outflow boundary. The calculated pressure and heat-transfer distributions 
are given in figures 7 and 8 where an excellent comparison with two other 
computational solutions is seen. The present analysis had no difficulty in 
handling this analytic body with its smooth variation in surface curvature. 
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The body-oriented coordinate system which has received very wide use 
because of its simplicity suffers several serious deficiencies, one of the 
more serious of which is that for bodies with a curvature discontinuity, such 
as spherically capped cones, there results a singularity in the governing 
equations. This problem can be circumvented by making the discontinuity 
point a special case (ref. 22) where special derivative relations are 
applied. However, the present numerically generated coordinate system does 
not suffer the same problem and this was amply demonstrated for inviscid flow 
(ref. 13) where solutions for a number of curvature discontinuous bodies were 
obtained for a wide range of conditions. Using the present analysis, the 
flow over a 45° spherically capped conical body was calculated using the 
condition of table 1. The computed heat-transfer distribution is plotted 
in figure 9 along with both experimental and computational results. Again, 
the comparisons are very good and the present analysis had no difficulty in 
handling the surface curvature discontinuity. 


A more rigorous demonstration of the capabilities of the numerical 
orthogonal coordinate generator can be seen on bodies with increasing nose 
bluntness, including those with reverse curvature. A family of bodies having 
these properties will be generated using the following cubic forebody 
generator: 

X = X + a«Y 2 + a_Y 3 
O 2 o 


where 


3 S l^l 

a 2 = V72 X 1 ' X o 


a 3 = - 2 a 2 V /3 V 


_ 1 

S 1 tan 0 C 
Xj = 0.5 

Y l = l. 

6 = 40° 

c 


The forebody joins smoothly to the conical afterbody with the surface 
angle specified by 0 c . Figure 10 gives the forebody shapes as a function 

of X Q and these five shapes represent the range of shapes that could be 

expected for planetary entry bodies experiencing severe ablation in the 
stagnation region (ref. 23). Using the conditions of table 1, but with 
a = 0.003, the flow fields about the five body shapes were computed. The 
pressure and heat- transfer distributions are given in figures 11 and 12 where 
the effect of nose blunting is readily apparent. Also, the representative 
effects of blunting on shock and sonic-line shape can be seen in figures 13-15 
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A representative converged coordinate mesh is given in figure 16 for the 

X =0.4 body. Note the highly compressed mesh near the body surface which 
o 

was necessary to resolve the boundary layer. All these results obtained 
with the present analysis are excellent with no noted undesirable flow-field/ 
coordinate system coupling effects. Part of the success of the present 
analysis is due to the smooth metric coefficients produced by the fourth- 
order differencing. Table 2 gives the metric coefficients at two locations 
for the X Q = 0.4 case and these results are typical of all the solutions 

obtained with the present analysis. 


CONCLUDING REMARKS 

The use of numerically generated orthogonal coordinates in compressible 
viscous-flow solutions to the Navier-Stokes equations is both practical 
and desirable. The technique has been demonstrated for a range of blunt - 
body shapes and the use of numerically evaluated metric coefficients has 
been successfully demonstrated. Although the present analysis had a moving 
boundary which necessitated recalculating the coordinate system at each step 
of the flow calculation, the technique should find easy application in 
static coordinate situations. 
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TABLE 1.- COMPUTATIONAL TEST CONDITIONS 


M m = 10.33 

CO 

Y = 1.1+ 

P = 100.77 N/m 2 
00 

T = U6.26° K 
00 

T = 330.6° K 
w 

^ = 0.03175 m 
N Re = 115600 
e = 2.5 x 10“ 6 
O = 0.001 


= 0.7 


TABLE 2.- METRIC COEFFICIENTS 

(a) Tangential Direction Variations 




TABLE 2.- METRIC COEFFICIENTS (CONCLUDED) 


(b) Normal Direction Variations 




0 ' 
.0333 
.0667 
.1000 
.1333 
.1667 
.2000 
.2333 
.2667 
.3000 
.3333 
.3267 
.Uooo 
.1*333 
.1:567 
.5000 
.5333 
.5667 
.6000 
.6333 

.6667 

.7000 

.7333 

.7667 

.8000 

.8333 

.8667 

.9000 

.9333 

.9667 

1.0000 


*1 

h 2 


.0021*9 

.00311 

1.6332 

.00389 

1.6331* 

.001*86 

1.6337 

.00608 

1.631*0 

.00759 

I.63UU 

'.0091*9 

1.631*9 

.01187 

1.6355 

.011*85 

1.6363 

.01856 

1.6373 

.02321 

1.6385 

.02902 

1.61*01 

.03629 

1.61*20 

.01*539 

1.61*1*1* 

.05679 

1.6471* 

.07105 

1.6512 

.08891 

1.6559 

.11129 

1.6617 

.13933 

1.6690 

.17UU9 

1.6781 

.21862 

I.689U 

.27398 

1.7035 

. 3l*3l*9 

1.7210 

.1*3075 

1.71*27 

. 5U021 

1.7696 

.67729 

1.8028 

.81*839 

1.81*33 

1.06071* 

1.8923 

1.32181 

1.9503 

1.63776 

2.0161* 

2.01539 



.1*529 

.U529 

.1*529 

.1*528 

.1*528 

.1*527 

.1*526 

.1*525 

.U52U 

.1*522 
.1*520 
.1*518 
.1*515 
.1*512 
.1*507 
.1*502 
.1*1*95 
. 1*U87 
.1*1*77 
. 1*1*65 
.1*1*53 
.1*1*30 
.1*1*22 
.1*1*07 
.1*395 
.1*38? 
.1*398 
.1*1*32 
.1*509 
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